data %>% glimpse
## Rows: 671
## Columns: 26
## $ birth <dbl> 81.511, 81.514, 81.552, 81.558, 81.593, 81.602, 81.610, NA, N…
## $ exit <dbl> 81.604, 81.539, 81.552, 81.667, 81.599, 81.771, 81.697, NA, N…
## $ hospstay <int> 34, 9, -2, 40, 2, 62, 32, NA, NA, 28, 38, NA, 62, 69, 1, 93, …
## $ lowph <dbl> NA, 7.250000, 7.059998, 7.250000, 6.969997, 7.189999, 7.32000…
## $ pltct <int> 100, 244, 114, 182, 54, NA, 282, NA, NA, 153, 229, NA, 182, 3…
## $ race <fct> white, white, black, black, black, white, black, NA, NA, blac…
## $ bwt <int> 1250, 1370, 620, 1480, 925, 940, 1255, 600, 700, 1350, 1310, …
## $ gest <dbl> 35.0, 32.0, 23.0, 32.0, 28.0, 28.0, 29.5, 26.0, 24.0, 34.0, 3…
## $ inout <fct> born at Duke, born at Duke, born at Duke, born at Duke, born …
## $ twn <int> 0, 0, 0, 0, 0, 0, 0, NA, NA, 0, 0, NA, 0, 0, 0, 0, 0, 1, 1, 0…
## $ lol <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ magsulf <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ meth <int> 0, 1, 0, 1, 0, 1, 1, NA, NA, 1, 0, NA, 0, 0, 1, 1, 1, 1, 1, 1…
## $ toc <int> 0, 0, 1, 0, 0, 0, 0, NA, NA, 0, 0, NA, 0, 0, 1, 1, 0, 0, 0, 1…
## $ delivery <fct> abdominal, abdominal, vaginal, vaginal, abdominal, abdominal,…
## $ apg1 <int> 8, 7, 1, 8, 5, 8, 9, NA, NA, 4, 6, NA, 6, 6, 2, 4, 8, 7, 1, 8…
## $ vent <int> 0, 0, 1, 0, 1, 1, 0, NA, NA, 0, 1, NA, 0, 0, 1, 1, 0, 0, 1, 1…
## $ pneumo <int> 0, 0, 0, 0, 1, 0, 0, NA, NA, 0, 0, NA, 1, 0, 1, 0, 0, 0, 1, 1…
## $ pda <int> 0, 0, 0, 0, 0, 0, 0, NA, NA, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0…
## $ cld <int> 0, 0, NA, 0, 0, 0, 0, NA, NA, 0, 0, NA, 1, 0, 0, 1, 0, 0, 0, …
## $ pvh <fct> NA, NA, NA, NA, definite, absent, NA, NA, absent, NA, NA, NA,…
## $ ivh <fct> NA, NA, NA, NA, definite, absent, NA, NA, absent, NA, NA, NA,…
## $ ipe <fct> NA, NA, NA, NA, NA, absent, NA, NA, absent, NA, NA, NA, absen…
## $ year <dbl> 81.51196, 81.51471, 81.55304, 81.55847, 81.59406, 81.60229, 8…
## $ sex <fct> female, female, female, male, female, female, female, NA, NA,…
## $ dead <int> 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0…
Сделайте копию датасета, в которой удалите колонки с количеством пропусков больше 100, а затем удалите все строки с пропусками.
cpdata <- data %>%
select(where(~sum(is.na(.))<=100)) %>%
drop_na()
cpdata %>% head
## birth exit hospstay lowph pltct race bwt gest inout twn
## 1 81.514 81.539 9 7.250000 244 white 1370 32.0 born at Duke 0
## 2 81.558 81.667 40 7.250000 182 black 1480 32.0 born at Duke 0
## 3 81.593 81.599 2 6.969997 54 black 925 28.0 born at Duke 0
## 4 81.610 81.697 32 7.320000 282 black 1255 29.5 born at Duke 0
## 5 81.624 81.700 28 7.160000 153 black 1350 34.0 born at Duke 0
## 6 81.626 81.730 38 7.039997 229 white 1310 32.0 born at Duke 0
## delivery apg1 vent pneumo pda cld year sex dead
## 1 abdominal 7 0 0 0 0 81.51471 female 0
## 2 vaginal 8 0 0 0 0 81.55847 male 0
## 3 abdominal 5 1 1 0 0 81.59406 female 1
## 4 vaginal 9 0 0 0 0 81.61053 female 0
## 5 abdominal 4 0 0 0 0 81.62421 female 0
## 6 vaginal 6 1 0 0 0 81.62695 male 0
Постройте графики плотности распределения для числовых переменных. Удалите выбросы, если таковые имеются. Преобразуйте категориальные переменные в факторы. Для любых двух числовых переменных раскрасьте график по переменной ‘inout’.
В оставшихся данных мы можем увидеть несколько категориальных переменных: twn (multiple gestation), vent (assisted ventilation used), pneumo (pneumothorax occurred), pda (patent ductus arteriosus detected), cld (on suppl. oxygen at 30 days), dead.
str(cpdata)
## 'data.frame': 531 obs. of 19 variables:
## $ birth : num 81.5 81.6 81.6 81.6 81.6 ...
## $ exit : num 81.5 81.7 81.6 81.7 81.7 ...
## $ hospstay: int 9 40 2 32 28 38 62 69 1 93 ...
## $ lowph : num 7.25 7.25 6.97 7.32 7.16 ...
## $ pltct : int 244 182 54 282 153 229 182 361 378 255 ...
## $ race : Factor w/ 4 levels "black","native American",..: 4 1 1 1 1 4 1 4 4 1 ...
## $ bwt : int 1370 1480 925 1255 1350 1310 1110 1180 970 770 ...
## $ gest : num 32 32 28 29.5 34 32 28 28 28 26 ...
## $ inout : Factor w/ 2 levels "born at Duke",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ twn : int 0 0 0 0 0 0 0 0 0 0 ...
## $ delivery: Factor w/ 2 levels "abdominal","vaginal": 1 2 1 2 1 2 2 1 2 2 ...
## $ apg1 : int 7 8 5 9 4 6 6 6 2 4 ...
## $ vent : int 0 0 1 0 0 1 0 0 1 1 ...
## $ pneumo : int 0 0 1 0 0 0 1 0 1 0 ...
## $ pda : int 0 0 0 0 0 0 0 0 0 0 ...
## $ cld : int 0 0 0 0 0 0 1 0 0 1 ...
## $ year : num 81.5 81.6 81.6 81.6 81.6 ...
## $ sex : Factor w/ 2 levels "female","male": 1 2 1 1 1 2 2 2 1 2 ...
## $ dead : int 0 0 1 0 0 0 0 0 1 0 ...
cpdata <- cpdata %>%
mutate_at(vars(twn, vent, pneumo, pda, cld, dead), as.factor)
Строим графики распределений:
distributions <- list()
nums <- cpdata %>% select(is.double | is.integer) %>% colnames()
idx <- 0
for (variable in nums) {
idx <- idx + 1
p <- ggplot(cpdata, aes_string(x=variable)) +
theme_minimal()
p <- if (idx < 3) {
p + geom_density(aes_string(fill='inout'), alpha=0.5) # ifelse почему то возвращает dataframe, а не ggplot объект. При этом при указании alpha в функции ggplot, параметр не распространяется на добавочные функции geom_density в рамках if else, странно
} else {
p + geom_density(fill='darkred', alpha=0.5)
}
distributions[[idx]] <- p
}
wrap_plots(distributions, axes='collect_y')
Попробуем удалить выбросы и посмотреть как изменятся распределения
remove_outliers <- function(data, column) {
Q1 <- quantile(data[[column]], 0.25)
Q3 <- quantile(data[[column]], 0.75)
IQR <- Q3 - Q1
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR
filtered_data <- data[data[[column]] >= lower_bound & data[[column]] <= upper_bound, ]
return(filtered_data)
}
cpdata_noout <- cpdata
for (variable in nums) {
cpdata_noout <- remove_outliers(cpdata_noout, variable)
}
distributions_noout <- list()
nums <- cpdata_noout %>% select(is.double | is.integer) %>% colnames()
idx <- 0
for (variable in nums) {
idx <- idx + 1
p <- ggplot(cpdata_noout, aes_string(x=variable)) +
theme_minimal()
p <- if (idx < 3) {
p + geom_density(aes_string(fill='inout'), alpha=0.5) # ifelse почему то возвращает dataframe, а не ggplot объект. При этом при указании alpha в функции ggplot, параметр не распространяется на добавочные функции geom_density в рамках if else, странно
} else {
p + geom_density(fill='darkred', alpha=0.5)
}
distributions_noout[[idx]] <- p
}
wrap_plots(distributions_noout, axes='collect_y')
В основном выбросы влияли на переменную hospstay, теперь она выглядит
намного лучше(отсутствуют отрицательные значений(скорее всего ошибка
ввода), график плотности стал менее вытянутым). Будем работать с данными
без выбросов.
Проведите тест на сравнение значений колонки ‘lowph’ между группами в переменной inout. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку ‘rstatix’. Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?
set.seed(123)
shapiro_test(cpdata_noout$lowph[cpdata_noout$inout == 'born at Duke'])
## # A tibble: 1 × 3
## variable statistic p.value
## <chr> <dbl> <dbl>
## 1 "cpdata_noout$lowph[cpdata_noout$inout == \"born at Duke\"]" 0.984 1.45e-4
shapiro_test(cpdata_noout$lowph[cpdata_noout$inout == "transported"])
## # A tibble: 1 × 3
## variable statistic p.value
## <chr> <dbl> <dbl>
## 1 "cpdata_noout$lowph[cpdata_noout$inout == \"transported\"]" 0.967 0.0562
if (shapiro_test(cpdata_noout$lowph[cpdata_noout$inout == 'born at Duke'])$p.value > 0.05 &&
shapiro_test(cpdata_noout$lowph[cpdata_noout$inout == "transported"])$p.value > 0.05) {
# Если данные нормально распределены, используем t-тест
test_result <- t_test(lowph ~ inout, data = cpdata_noout)
print('Использован Т-тест')
} else {
# Если данные не нормально распределены, используем U-тест Манна-Уитни
test_result <- wilcox_test(lowph ~ inout, data = cpdata_noout)
print('Использован тест Манна-Уитни')
}
## [1] "Использован тест Манна-Уитни"
# Выводим результаты теста
print(test_result)
## # A tibble: 1 × 7
## .y. group1 group2 n1 n2 statistic p
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl>
## 1 lowph born at Duke transported 413 73 20657 0.000000446
ggplot(cpdata_noout, aes(x = inout, y = lowph)) +
geom_boxplot() +
stat_compare_means(method = "wilcox.test") +
labs(title = "Сравнение значений lowph между группами",
x = "Группы",
y = "Значения lowph")
Интерпретация: тк среднее значение low pH статистически значимо ниже у транспортированных пациентов и его снижение положительно ассоциированно с риском смерти, то можно предположить, что транспортировка (условия транспортировки и/или время затраченное на нее) ухудшают состояние пациента и снижают его шансы на выживание.
Сделайте новый датафрейм, в котором оставьте только континуальные или ранговые данные, кроме ‘birth’, ‘year’ и ‘exit’. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.
ftask <- cpdata_noout %>% select(-c(birth, year, exit)) %>% select(is.double|is.factor)
cormat <- ftask %>% mutate(across(where(is.factor), as.numeric)) %>% cor()
corrplot(cormat)
ggcorrplot(cormat)
Постройте иерархическую кластеризацию на этом датафрейме.
ftask <- cpdata_noout %>% select(-c(birth, year, exit)) %>% select(is.double|is.factor) %>% mutate(across(where(is.factor), as.numeric))
distmat <- dist(t(scale(ftask)), method = "euclidean")
hc <- hclust(distmat, method = "complete")
plot(hc, main = "Hierarchical Clustering Dendrogram", xlab = "Variables", sub = "", cex = 0.8)
pheatmap(
scale(ftask),
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "complete",
main = "Heatmap with Hierarchical Clustering"
)
Интерпретация:
Группа умерших людей подсвечена красным в колонке dead. По хитмапу можно сказать, что pda, pneumo, vent имеют хорошие перекрытия с dead (и являются частью одной клады), при этом lowph тоже имеют хорошую отрицательную корреляцию с исходом, эти переменные можно было бы попробовать использовать в качестве предикторов.
Проведите PCA анализ на этих данных. Проинтерпретируйте результат. Нужно ли применять шкалирование для этих данных перед проведением PCA?
Переменные в нашем датасете имеют разные диапозоны значений, так что шкалирование необходимо, чтобы исключить доминирование “больших” переменных над малыми
df_scaled <- scale(ftask)
pca_result <- prcomp(df_scaled, center = TRUE, scale. = TRUE)
summary(pca_result)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.7509 1.1969 1.05804 1.04516 0.97608 0.93719 0.90775
## Proportion of Variance 0.2555 0.1194 0.09329 0.09103 0.07939 0.07319 0.06867
## Cumulative Proportion 0.2555 0.3749 0.46815 0.55918 0.63857 0.71176 0.78043
## PC8 PC9 PC10 PC11 PC12
## Standard deviation 0.83482 0.76392 0.74432 0.66321 0.60038
## Proportion of Variance 0.05808 0.04863 0.04617 0.03665 0.03004
## Cumulative Proportion 0.83851 0.88714 0.93331 0.96996 1.00000
Видно, что для объяснения 75 процентов дисперсии необходимо взять целых 7(из 12) компонент, значит, что в данных нету переменных, которые бы однозначно их описывали.
# Выполнение PCA
pca_result <- prcomp(df_scaled, center = TRUE, scale. = TRUE)
# Построение biplot
p <- fviz_pca_biplot(
pca_result,
text = paste0('Patient ID: ', rownames(ftask)),
geom.ind = "point",
palette = c("blue", "red"),
repel = TRUE
) + labs(title = "PCA Biplot") +
geom_point(aes(col= as.factor(ftask$dead)), text = paste0('Patient ID: ', rownames(ftask))) +
scale_color_discrete('Patient', labels = c('1' = 'Alive', '2' = 'Dead'))
p
Переведите последний график в ‘plotly’. При наведении на точку нужно, чтобы отображалось id пациента.
plot_ly(
x = pca_result$x[,1],
y = pca_result$x[,2],
text = pca_result$x %>% rownames(),
hovertemplate = paste('<br><b>Patient ID</b>: <b>%{text}</b>'),
hoverinfo = 'text',
type = 'scatter',
mode = 'markers',
color = fct_recode(as.factor(ftask$dead), 'Alive' = '1', 'Dead' = '2')
) %>%
layout(title = "PCA Biplot",
xaxis = list(title = "Dim2(11%)"),
yaxis = list(title = "Dim1(23.6%)"),
legend = list(title = list(text = "Patient Status")))